function [pathwayStatistics,listOfPoolMetabolites,zeroSolution, ...
    listOfPathways,listOfAllWeights, zeroSolutionBiomass, modelToReturn]=...
        get_pathway_statistics(pathway, model, minimumBiomassFlux, shouldModelBeReturned,...
        defines, thermodinamicsValues)

     if strcmp(defines.isThermodinamicUsed, 'true')
        thermodinamicsValues = evalin('base', 'thermodinamicsValues'); 
     end
    
    currentModel = model;
    addedPoolMetabolites = [];
    allMetabolitesInPathway = [];
    addedMetabolites = [];
    existingMetabolites = [];
    listOfPoolMetabolites = [];
    sourceTransportReaction ='';
    zeroSolution = false;
    zeroSolutionBiomass = false;
    source_metabolite = pathway(1,:);
    listOfPathways = {};
    modelToReturn = '';
    
    pathwayStatistics = struct( ...
                            'MET_COUNT',0, ...
                            'EXIST_MET_COUNT',0, ...
                            'ADDED_MET_COUNT',0, ...
                            'REAC_COUNT',0, ...
                            'EXISTING_REACT_COUNT',0, ...
                            'REACTIONS_WITH_KNOWN_THERMODINAMICS',0, ...
                            'UNBALANCED_REACT_COUNT',0, ...
                            'FBA_RESULT', 0, ...   
                            'FBA_BIOMASS_RESULT', 0 );

    for y = 1:length(pathway(:,1))
       if regexp(pathway(y,:), 'R') == 1
            listOfAnalyzedMetabolites = [];
            pathwayStatistics.REAC_COUNT = pathwayStatistics.REAC_COUNT + 1;
            reactionStructure = struct;
            equation = ''; 
            reactionStructure = get_reaction_information(pathway(y,:));

            for j = 1:length(reactionStructure.PATHWAY)
                if length(listOfPathways) > 0
                    [isReactionMemeberOfPathways, index] = ismember(reactionStructure.PATHWAY{j}, listOfPathways(:,1));
                else
                    isReactionMemeberOfPathways = false;
                end

                if isReactionMemeberOfPathways == 0
                   listOfPathways = [listOfPathways; {reactionStructure.PATHWAY{j}}, 1];
                else
                   listOfPathways(index,2) = {listOfPathways{index,2} + 1};
                end
            end

            equation = getfield(reactionStructure, 'EQUATION');              
            metabolites = regexp(equation, 'C.....', 'match');

            metabolitesCountInReaction = 0;
            metaboliteCountAddedToModel = 0;
            metabolitesCountExistingInModel = 0;

            for n = 1:size(metabolites,2)
                metaboliteStructure = struct;

               [metaboliteAbbreviation,isFoundInModel] =  ...
                   get_metabolite_abbreviation_from_model(currentModel, metabolites{n}, defines.cystoleCompartment);

               temporaryMetabolite = metabolites{n};
                if ~isFoundInModel
                    addedMetabolites = [addedMetabolites; metabolites{n}];
                end
                if isempty(strmatch(metabolites{n}, pathway))
                    allMetabolitesInPathway = [allMetabolitesInPathway;metabolites{n}];
                end
                if length(metaboliteAbbreviation) > 12
                    if regexp(equation, metaboliteAbbreviation(1:12))
                        temporaryMetabolite = metaboliteAbbreviation;
                    end
                end
                equation = strrep(equation, temporaryMetabolite, metaboliteAbbreviation);
                metabolitesCountInReaction = metabolitesCountInReaction + 1;
                if isFoundInModel == 1
                    if isempty(~strmatch(metabolites{n},existingMetabolites))
                       if isempty(regexp(metaboliteAbbreviation, 'added', 'once'))
                           metabolitesCountExistingInModel = metabolitesCountExistingInModel + 1;           
                           existingMetabolites = [existingMetabolites; metabolites{n}];
                       end
                    end
                else
                    metaboliteCountAddedToModel = metaboliteCountAddedToModel + 1;
                end
                metaboliteStructure = get_compound_information(metabolites{n});
                listOfAnalyzedMetabolites = [listOfAnalyzedMetabolites; {metaboliteAbbreviation, metaboliteStructure}];
            end
            
            reactionStructure = setfield(reactionStructure,'EQUATION_FROM_MODEL', equation);
            reactionStructure = setfield(reactionStructure,'metabolitesCountInReaction', metabolitesCountInReaction);
            reactionStructure = setfield(reactionStructure,'metabolitesCountExistingInModel', metabolitesCountExistingInModel);
            reactionStructure = setfield(reactionStructure,'metaboliteCountAddedToModel', metaboliteCountAddedToModel);

            name = getfield(reactionStructure, 'NAME');
            reactionStructure.NAME = regexprep(name,' ','_');
            pathwayStatistics.MET_COUNT = pathwayStatistics.MET_COUNT + metabolitesCountInReaction;
            pathwayStatistics.EXIST_MET_COUNT = pathwayStatistics.EXIST_MET_COUNT + metabolitesCountExistingInModel;
            pathwayStatistics.ADDED_MET_COUNT = pathwayStatistics.ADDED_MET_COUNT + metaboliteCountAddedToModel;

            currentReactionEquation = equation;
            reactionsNumberWithThermodinamics = 0;
            
            if strcmp(defines.isThermodinamicUsed,'true')
                isThermodinamicsValueFound = is_thermodinamics_data_found(pathway(y,:),thermodinamicsValues);
                reactionsNumberWithThermodinamics = reactionsNumberWithThermodinamics + isThermodinamicsValueFound;
            end
            
            reactionName = reactionStructure.NAME;
            [currentModel, iFExist] = addReaction(currentModel, {pathway(y,:),reactionName}, currentReactionEquation); 
                
            if iFExist > 0
                pathwayStatistics.EXISTING_REACT_COUNT = pathwayStatistics.EXISTING_REACT_COUNT + 1;
            end

            for z = 1:length(listOfAnalyzedMetabolites)
                if regexp(listOfAnalyzedMetabolites{z,1}, 'added_')
                    if regexp(listOfAnalyzedMetabolites{z,1}, source_metabolite)
                        sourceTransportReaction = strcat('R_', source_metabolite, '_transport');
                        if ~ismember(sourceTransportReaction, currentModel.rxns)                       
                           exctracellularSourceMetabolite = strcat('added_',source_metabolite, defines.extracellularCompartment);                   
                           transportReactionEquation = sprintf('%s %s %s',listOfAnalyzedMetabolites{z,1}, '->', exctracellularSourceMetabolite);
                           currentModel = addReaction(currentModel, sourceTransportReaction, transportReactionEquation); 

                           sourceMetaboliteExchangeReaction = strcat('EX_',source_metabolite,'(e)');
                           exchangeReactionEquation = sprintf('%s %s %s',exctracellularSourceMetabolite, '<=>');
                           currentModel = addReaction(currentModel, sourceMetaboliteExchangeReaction, exchangeReactionEquation); 

                           sourceMetaboliteIndex = find(strcmp(listOfAnalyzedMetabolites{z,1}, currentModel.mets),1);
                           sourceExtracellularMetaboliteIndex = find(strcmp(exctracellularSourceMetabolite, currentModel.mets),1);

                           currentMetaboliteStructure = listOfAnalyzedMetabolites(z, 2);

                           currentModel.metFormulas{sourceMetaboliteIndex} = currentMetaboliteStructure{1}.FORMULA;
                           currentModel.metFormulas{sourceExtracellularMetaboliteIndex} = currentMetaboliteStructure{1}.FORMULA;

                           keggID = regexp(currentMetaboliteStructure{1}.ENTRY, 'C.....', 'match');
                           currentModel.metKEGGID{sourceMetaboliteIndex} = keggID{1};
                           currentModel.metKEGGID{sourceExtracellularMetaboliteIndex} = keggID{1};
                        end
                    else
                        pool_metabolite_transport_rxn = strcat('R_', listOfAnalyzedMetabolites{z,1},'_transport');
                        if ~ismember(pool_metabolite_transport_rxn, currentModel.rxns)  
                            isMetaboliteFoundInPathway = false;
                            for x = 1:length(pathway(:,1))
                                if regexp( listOfAnalyzedMetabolites{z,1}, pathway(x,:))>0
                                   isMetaboliteFoundInPathway = true;
                                end
                            end
                            currentMetaboliteStructure = listOfAnalyzedMetabolites(z, 2);
                            if ~isMetaboliteFoundInPathway
                                extracellulatrPoolMetabolite = strcat(strrep(listOfAnalyzedMetabolites{z,1}, defines.cystoleCompartment, defines.extracellularCompartment));
                                transportReactionEquation = sprintf('%s %s %s',listOfAnalyzedMetabolites{z,1}, '<=>', extracellulatrPoolMetabolite);
                                currentModel = addReaction(currentModel, pool_metabolite_transport_rxn, transportReactionEquation); 

                                poolMetaboliteExchangeReactionEquation = strcat('EX_',strrep(listOfAnalyzedMetabolites{z,1},'[c]','(e)')); 
                                exchangeReactionEquation = sprintf('%s %s',extracellulatrPoolMetabolite, '<=>' );
                                currentModel = addReaction(currentModel, poolMetaboliteExchangeReactionEquation,exchangeReactionEquation); 

                                index_of_pool_extracellular_metaboilite = find(strcmp(extracellulatrPoolMetabolite, currentModel.mets),1);
                             
                                keggID = regexp(currentMetaboliteStructure{1}.ENTRY, 'C.....', 'match');
                                currentModel.metFormulas{index_of_pool_extracellular_metaboilite} = currentMetaboliteStructure{1}.FORMULA;
                                
                                if isempty(keggID)
                                    
                                else
                                currentModel.metKEGGID{index_of_pool_extracellular_metaboilite} = keggID{1};
                                end
                                pool_met =  regexprep(regexprep(listOfAnalyzedMetabolites{z,1}, '[[c]]',''),'added_','');
                                
                                listOfPoolMetabolites = [listOfPoolMetabolites; pool_met];
                                if ~regexp('*', listOfAnalyzedMetabolites{z,2}.ENTRY)
                                    addedPoolMetabolites = [addedPoolMetabolites; listOfAnalyzedMetabolites{z,2}.ENTRY];
                                end
                            end
                            
                            poolMetaboliteIndex = find(strcmp(listOfAnalyzedMetabolites{z,1}, currentModel.mets),1);
                            currentModel.metFormulas{poolMetaboliteIndex} = currentMetaboliteStructure{1}.FORMULA;
                            keggID = regexp(currentMetaboliteStructure{1}.ENTRY, 'C.....', 'match');
                            if ~isempty(keggID)
                            currentModel.metKEGGID{poolMetaboliteIndex} = keggID{1};
                            end
                        end
                    end
                end
            end

           if ~iFExist
                addedReactionIndex = find(strcmp(pathway(y,:), currentModel.rxns),1);
                [unbalancedRxns] = checkMassChargeBalance(currentModel,addedReactionIndex);

                if unbalancedRxns > 0
                    pathwayStatistics.UNBALANCED_REACT_COUNT = pathwayStatistics.UNBALANCED_REACT_COUNT + 1;
                end       
           end
        elseif regexp(pathway(y,:), 'C') == 1
            allMetabolitesInPathway = [allMetabolitesInPathway; pathway(y,:)];
        end
    end     

    solution = optimizeCbModel(currentModel);
    pathwayStatistics.FBA_BIOMASS_RESULT = solution.f;
    
    if round(solution.f) <= 0
        zeroSolutionBiomass = true;
    end

    objectiveReaction = currentModel.rxns(find(currentModel.c,1));
    currentModel = changeRxnBounds(currentModel, objectiveReaction, minimumBiomassFlux,'l');
    currentModel = changeObjective(currentModel, sourceTransportReaction);
    solution = optimizeCbModel(currentModel);
    pathwayStatistics.FBA_RESULT = solution.f;

    if round(solution.f) <= 0
        zeroSolution = true;
    end
    
    pathwayStatistics.MET_COUNT = pathwayStatistics.MET_COUNT - (((length(pathway(:,1)) - 1) / 2) - 1);
    pathwayStatistics.REACTIONS_WITH_KNOWN_THERMODINAMICS = reactionsNumberWithThermodinamics;
    reactions_weight = calculate_reactions_weight( ...
        pathwayStatistics.REAC_COUNT, ...
        pathwayStatistics.EXISTING_REACT_COUNT, ...
        pathwayStatistics.REAC_COUNT - pathwayStatistics.UNBALANCED_REACT_COUNT, ...
        pathwayStatistics.REACTIONS_WITH_KNOWN_THERMODINAMICS);
    all_metabolites_weight=0;
    
    for i = 1:length(allMetabolitesInPathway(:,1))
        weight = calculate_metabolites_occurence_weight(allMetabolitesInPathway(i,:), currentModel);
        all_metabolites_weight = all_metabolites_weight+weight;
    end
    
    all_metabolites_weight = all_metabolites_weight / length(allMetabolitesInPathway(:,1));
    existing_metabolites_weight = calculate_existing_metabolites_weight(pathwayStatistics.EXIST_MET_COUNT, pathwayStatistics.MET_COUNT);

    if isempty(listOfPoolMetabolites)
        count_of_pool_metabolites = 0;
    else
        count_of_pool_metabolites =  length(listOfPoolMetabolites(:,1));
    end

    added_metabolites_weight = calculate_added_metabolites_weight( pathwayStatistics.REAC_COUNT,count_of_pool_metabolites, ...
        pathwayStatistics.ADDED_MET_COUNT,pathwayStatistics.MET_COUNT);

    listOfAllWeights = [round(pathwayStatistics.FBA_RESULT*100)/100, ...
                        round(pathwayStatistics.FBA_BIOMASS_RESULT*100)/100, ...
                        round(reactions_weight*100)/100, ...
                        round(all_metabolites_weight*100)/100 + ... 
                        round(existing_metabolites_weight*100)/100 + ...
                        round(added_metabolites_weight*100)/100];
                    
    if shouldModelBeReturned
        modelToReturn = currentModel;
    end
end



